##############################################################################
##### Robustness Test Legitimation for QCA on Authoritarian Persistence  #####
##############################################################################


rm(list = ls())

# set working directory:
setwd()

## NOTE: For correct replications, please load the packages QCA 2.2 and SetMethods 2.0

# load packages
library(QCA); library(QCAGUI); library(SetMethods); library(lattice);
library(arm); library(plyr); library(car); library(stringr); library(xtable)


# load data from the following sources: -----------------------------------------------#
# -------------------------------------------------------------------------------------#


# 1.) for outcome "authoritarian persistence": 
# Geddes, Barbara, Erika Frantz, and Joseph Wright. 2014. 
# "Autocratic Breakdown and Regime Transitions: A New Data Set." 
# Perspectives on Politics 12(2): 313-31. 

# 2.) for both forms of repression:
# Cingranelli, David L, and David L Richards. 2013. 
# "The CIRI Human Rights Dataset, Version 2012.04.10." 
# http://www.humanrightsdata.com/ (May 21, 2016).

# 3.) for cooptation as compensating vulnerability:
# Schmotz, Alexander. 2015. "Vulnerability and Compensation: 
# Constructing an Index of Co-Optation in Autocratic Regimes." 
# European Political Science 14(4): 439-57. 

# 4.) for cooptation as simulating pluralism:
# Beck, Thorsten et al. 2001. 
# "New Tools in Comparative Political Economy: 
# The Database of Political Institutions." 
# The World Bank Economic Review 15(1): 165-76.

# 5.) for specific and diffuse support of legitimation:
# Grauvogel, Julia, and Christian von Soest.  
# "Unpacking Authoritarian Legitimation Strategies: 
# Findings from a New Explorative Expert Survey." 
# Forthcoming in Zeitschrift f?r Vergleichende Politikwissenschaft.


# prepare and adjust raw data as explained in paper -------------------------------------#
# (diverging regime spans, adjust period of observation, apply weighted average) --------#


# calibrate conditions with alternative anchors for legitimation -----------------#
# --------------------------------------------------------------------------------#

# Set calibration (using loop function): 
# STEP 1:
# store all of your conditions in an object 

conds <- c("repp","repc","coopv","fomyth", "ideo", "perso", "inteng", "proc", "perf")


# STEP 2:
# create matrix that stores the qualitative anchors for each condition (cf. paper)
# the order of conditions must be the one in which they appear in object conds

conds_anchors <- matrix(c(-1, 1.2, 5,
                          6, 4.2, 1,
                          9, 5.4, 2,
                          1, 2.99, 5,
                          1, 2.77, 4.25,
                          1, 2.84, 4.25,
                          1, 1.9, 4,
                          1, 3.4, 5,
                          1, 3.1, 4.75),
                        ncol=3, byrow = TRUE)

# NOTE: This robustness test uses the anchors as suggeste by Grauvogel and von Soest
# to see how this affects the QCA results


# STEP 3:
# now use those qualitative anchors stored 
# in object conds_anchors in the calibration function
# note that we store the newly calibrated sets in our data set, which happens to have 15
# columns, so we start putting the new sets at 15+i, 

for (i in 1:length(conds)) {
  data[ ,15+i] <- calibrate(data[ ,conds[i]], type = "fuzzy", logistic = TRUE, idm = 0.95,
  thresholds = c("e"= conds_anchors[i,1],"c"= conds_anchors[i,2],"i"= conds_anchors[i,3]))
}


# STEP 4
# name the new columns just added to the data, using their original name with "high" 
names(data)[c(16:24)] <- paste0("high", conds)


# calibrate crisp set on coops (cf. paper)
# set anchors in new vector (coop1)
coop <- subset(data, data$coops > 4)
coop1 <- subset(coop, coop$coops < 7)

# write function for matching the country names in data and coop1 
# and code the matchings as 1 (highcoops) and 0 (no highcoops) in new column

# create n (lenght of rows)
n <- nrow(data) 
for (i in 1:n) 
{
  str <- data$cowcode[i]
  index <- pmatch(str, coop1$cowcode)
  if (is.na(index)){
    data$highcoops[i] <- 0
  }
  else{
    data$highcoops[i] <- 1
  }
}

# recode Iran, based on substantive arguments (cf. paper)
# it is considered to be out of the set "coops" despite its DPI scores
# (the condition coops is in column 25)
data["IRN", 25] <- 0


# STEP 5
# Create histograms and plots using a loop function

# create object with all sets for histograms (including highcoops)
setsclo <- c("highrepp", "highrepc","highcoops","highcoopv","highfomyth", "highideo", "highperso", 
             "highinteng", "highproc", "highperf")
setsclo

# histograms for all calibrated sets
par(mfrow=c(3,4))
for (i in setsclo){
  hist(data[, i],
       xlab=i, main = paste(""))
}

# create object with all fuzzy sets for plots (excluding highcoops)
sets <- c("highrepp", "highrepc","highcoopv","highfomyth", "highideo", "highperso", 
          "highinteng", "highproc", "highperf")
sets

# plots of sets againts base variables
par(mfrow=c(3,4))
for (i in 1:length(conds)) {
  plot(data[, conds[i]], data[, sets[i]],
       xlab=conds[i], ylab = sets[i], 
       main = paste(" "))
  abline(h=0.5, col="black")
  abline(v=conds_anchors[i,2]) 
}

# xy plots each set and outcome
par(mfrow=c(3,4))
for (i in sets) {
  xy.plot(data[,i], data$pers,
          main = i,
          xlab = i,
          ylab = 'persistent')
}

# for alternative calibrations see other scripts on robustness tests


# apply hexagon ----------------------------------------------------------------#
# ------------------------------------------------------------------------------#

# to prevent confusion, erase uncalibrated columns:
data[ ,(6:15)] <- NULL

# diffuse support of legitimation, fomyth, ideo and perso combined by OR
data$legd <- pmax(data$highfomyth, data$highideo, data$highperso)

# specific support of legitimation, inteng, proc,perf combined by OR
data$legs <- pmax(data$highinteng, data$highproc, data$highperf)

# exclude not merged leg data
data[ ,c(9:14)] <- NULL

# rename other colums:
data <- rename(data, c("highrepp" = "repp", 
                       "highrepc" = "repc",
                       "highcoopv" = "coopv",
                       "highcoops" = "coops"))



# checking the effects of alternative anchors on the QCA results -------------------------#
# Outcome: persistent regimes ------------------------------------------------------------#


# Analysis of necessity ------------------------------------------------------------------#


QCAfit(data[, 6:11], data$pers, necessity = TRUE, names(data[, 6:11]))
QCAfit(1-data[, 6:11], data$pers, necessity = TRUE, paste("not", names(data[, 6:11])))

SUIN <- superSubset(data, outcome = "pers",
                    conditions = c("repp", "repc", "legd", "legs","coopv", "coops"),
                    incl.cut = 0.9, cov.cut = 0.7,
                    necessity = TRUE)
SUIN


# Analysis of sufficiency ----------------------------------------------------------------#


# truth Table
TT_y <- truthTable(data, outcome = "pers",
                   conditions = c("repp", "repc", "legd", "legs","coopv", "coops"),
                   incl.cut1 = 0.77,
                   n.cut = 1,
                   show.cases = TRUE,
                   details = TRUE,
                   complete = FALSE,
                   sort.by = c("OUT", "incl", "n"))

TT_y




# conservative solution
sol_y_c <- eqmcc(data, outcome = "pers",
                 conditions = c("repp", "repc", "legd", "legs","coopv", "coops"),
                 incl.cut1 = 0.77,
                 n.cut = 1,
                 show.cases = TRUE,
                 details = TRUE)
sol_y_c


# most parsimonious solution
sol_y_p <- eqmcc(data, outcome = "pers",
                 conditions = c("repp", "repc", "legd", "legs","coopv", "coops"),
                 incl.cut1 = 0.77,
                 n.cut = 1,
                 show.cases = TRUE,
                 details = TRUE,
                 row.dom = TRUE,
                 include = "?")
sol_y_p

# simplifying assumptions
sol_y_p$SA

# intermediate solution
sol_y_i <- eqmcc(data, outcome = "pers",
                 conditions = c("repp", "repc", "legd", "legs","coopv", "coops"),
                 incl.cut1 = 0.77,
                 n.cut = 1,
                 show.cases = TRUE,
                 details = TRUE,
                 include = "?",
                 dir.exp = c("-",1,1,1,1,1))
sol_y_i

# dir. ex.: all conds are expected to be present in all persistent regimes beside repp due 
# to its backfiring effects (cf. paper)

# easy counterfactuals
sol_y_i$i.sol$C1P1$EC


# summary of robustness test 3:
# Grauvogel and von Soest mostly took the mean values (with small deviations) 
# for determining the cross over point,
# this can result in flawed findings (Schneider and Wagemann, 2012, 33); 
# furthermore, the plots illustrate
# that the distribution of all legitimation fuzzy scores is skewed;
# concerning the QCA results, these anchors cause model ambiguity and less coverage
# compared to my solution formulas

